from pylab import *
from scipy.ndimage import measurements
import statsmodels.api as sm

L = 150
r = rand(L,L)
p = 0.4
z = r < p

imshow(z, origin='lower', interpolation='nearest')
colorbar()  #增加颜色类标
title("Matrix")
plt.show()  #很关键


lw, num = measurements.label(z)
imshow(lw, origin='lower', interpolation='nearest')
colorbar()
title("Labeled clusters")
plt.show()  #很关键



b = arange(lw.max() + 1)
shuffle(b)
shuffledLw = b[lw]
imshow(shuffledLw, origin='lower', interpolation='nearest')
colorbar()
title("Labeled clusters")
plt.show()  #很关键



area = measurements.sum(z, lw, index=arange(lw.max() + 1))
areaImg = area[lw]
im3 = imshow(areaImg, origin='lower', interpolation='nearest')
colorbar()
title("Clusters by area")
plt.show()  #很关键


im3 = imshow(areaImg, origin='lower', interpolation='nearest')
colorbar()
title("Clusters by area")
sliced = measurements.find_objects(areaImg == areaImg.max())
if(len(sliced) > 0):
    sliceX = sliced[0][1]
    sliceY = sliced[0][0]
    plotxlim=im3.axes.get_xlim()
    plotylim=im3.axes.get_ylim()
    plot([sliceX.start, sliceX.start, sliceX.stop, sliceX.stop, sliceX.start], \
    [sliceY.start, sliceY.stop, sliceY.stop, sliceY.start, sliceY.start], \
    color="red")
    xlim(plotxlim)
    ylim(plotylim)



def RankOrderPlot(data):
    d=array(data)
    d = d[d>0]
    t=array(sorted(d,key=lambda x:-x))
    r=array(range(1,len(d)+1))
    y = log(t)
    x = log(r)
    X = sm.add_constant(x, prepend=True)
    res = sm.OLS(y,X).fit()
    C,beta = res.params
    plt.plot(r,t,"o",color="b")
    plt.plot(r,exp(C)*r**beta,"r-")
    plt.xscale('log')
    plt.yscale('log')
    plt.text(max(r)/2,max(t)/100,"beta = " + str(round(beta,2)))
    print (beta)
RankOrderPlot(area)
plt.show()  #很关键